# Who Cares if You Vote? Partisan agreement and injunctive norms of voting
# Edward Fieldhouse (1), David Cutts (2), Jack Bailey (1)
# (1) The University of Manchester, (2) The University of Birmingham


# 1. Housekeeping ---------------------------------------------------------

# The following code requires that you have first opened the associated
# project file. If not, click the button on the top right and open it
# ("Open Project...").

# Note also that the Bayesian methods we use below rely on having installed
# the probabilistic programming language Stan. If you have yet to install
# Stan, please see https://mc-stan.org/users/interfaces/.

# Load packages

library(tidyverse)
library(tidybayes)
library(magrittr)
library(ggpubr)
library(jbmisc)   # devtools::install_github("jackobailey/jbmisc)
library(haven)
library(brms)
library(here)


# Load data from the British Election Study Internet Panel. Note: If you
# plan to use this data, please make sure that you download the most recent
# version available at https://www.britishelectionstudy.com.

bes <- read_sav(here("_data", "BES2017_W14_Panel_v0.4.sav.zip"))



# 2. Transform data (respondent-level) ------------------------------------

# Select variables

bes <-
  bes %>%
  select(
    # Respondent-level covariates
    id,
    wave12,
    age = Age,
    female = gender,
    edu = edlevelW12,
    edu2 = educationW12,
    married = marital,
    pid = partyIdW12,
    sq = partyIdSqueezeW12,
    str = partyIdStrengthW12,
    
    # Discussant-level covariates
    rel_1 = relationshipName1W12,
    rel_2 = relationshipName2W12,
    rel_3 = relationshipName3W12,
    app_1 = discussantApprovalVoteName1W12,
    app_2 = discussantApprovalVoteName2W12,
    app_3 = discussantApprovalVoteName3W12,
    vote_1 = discussantVoteName1W12,
    vote_2 = discussantVoteName2W12,
    vote_3 = discussantVoteName3W12) %>%
  filter(wave12 == 1)


# Z-standardise age

bes$age_z <-
  bes$age %>% 
  rescale()


# Convert female to factor with male as reference

bes$female <-
  bes$female %>%
  as_factor()


# Add "other" and "unknown" categories to education variable

bes$edu <-
  bes$edu %>%
  as_factor() %>%
  factor(levels = c("No qualifications",
                    "Below GCSE",
                    "GCSE",
                    "A-level",
                    "Undergraduate",
                    "Postgrad",
                    "Other",
                    "Unknown"))


# Populate "other" and "unknown" education categories from educationW1_6 variable

bes$edu[bes$edu2 == 18] <- "Other"
bes$edu[is.na(bes$edu) == T] <- "Unknown"


# Marital status

bes$married <-
  bes$married %>%
  as_factor() %>%
  as_dummy(c("Married", "Living as married", "Civil Partnership")) %>%
  factor(labels = c("Other",
                    "Cohabiting"))


# PID: Combine PID and PID Squeeze

bes$pid[bes$pid %in% c(10, 9999)] <- bes$sq[bes$pid %in% c(10, 9999)]


# Party ID

bes$pid <-
  bes$pid %>%
  as_factor() %>%
  factor(labels = c("Conservative",
                    "Labour",
                    "Liberal Democrat",
                    "SNP",
                    "Plaid Cymru",
                    "UKIP",
                    "Green Party",
                    "Other",
                    "None",
                    "Don't know"))

bes$sq <- NULL


# PID Strength

bes$str <- 
  bes$str %>% 
  as_factor() %>% 
  mark_na("Don't know") %>% 
  fct_explicit_na("None/Don't know") %>% 
  fct_relevel("None/Don't know")


# Convert data from respondent-level to dyad-level

bes <- 
  bes %>% 
  data.frame() %>% 
  reshape(
    varying = names(bes)[10:18],
    sep = "_",
    ids = "id",
    direction = "long",
    timevar = "disc"
  )



# 3. Transform data (referent-level) --------------------------------------

# Dyadic relationship, reference = "Other"

bes$rel <- 
  bes$rel %>% 
  as_factor() %>%
  fct_relevel("Other")


# Dyadic approval

bes$app <- 
  bes$app %>% 
  as_dummy(1) %>% 
  factor(labels = c("Other", "Yes"))


# Dyadic partisanship

bes$vote <- 
  bes$vote %>% 
  as_factor() %>%
  factor(labels = c("Labour",
                    "Conservative",
                    "Liberal Democrat",
                    "SNP",
                    "Plaid Cymru",
                    "Green Party",
                    "UKIP",
                    "Other",
                    "None",
                    "Don't know"))

opt <- c("None", "Don't know")

bes <- 
  bes %>% 
  mutate(
    prt =
      case_when(
        pid == vote & !(pid %in% opt) & !(vote %in% opt) & pid != "Other" & vote != "Other" ~ "Shared Partisanship",
        pid != vote & !(pid %in% opt) & !(vote %in% opt) ~ "Opposing Partisanship",
        pid == "Other" & vote == "Other" ~ "Opposing Partisanship",
        pid %in% opt & !(vote %in% opt) ~ "Respondent Non-partisan",
        vote %in% opt & !(pid %in% opt) ~ "Discussant Non-partisan",
        pid %in% opt & vote %in% opt ~ "Shared Non-partisan"
      ) %>% 
      factor(levels = c("Shared Non-partisan",
                        "Shared Partisanship",
                        "Opposing Partisanship",
                        "Respondent Non-partisan",
                        "Discussant Non-partisan"))
  )


# Remove NAs listwise

bes <-
  bes %>%
  na.omit()


# Compute number of discussants by ID

bes <- 
  bes %>% 
  add_count(id,
            name = "n_disc")



# 4. Fit model ------------------------------------------------------------

# Fit multilevel approval model with brms

m1 <- brm(formula = 
            app ~ 1 + rel + prt + age_z + female + edu + married + str + n_disc + (1 | id),
          family = bernoulli(link = "logit"),
          prior = 
            prior(normal(0, 1.5), class = "Intercept") +
            prior(normal(0, 0.5), class = "b") +
            prior(exponential(2), class = "sd"),
          data = bes,
          iter = 2000,
          chains = 4,
          cores = 4,
          refresh = 5,
          control = list(adapt_delta = .99,
                         max_treedepth = 15),
          file = here("_output", "app1_w12"))


# Inspect model

plot(m1, ask = FALSE)
summary(m1)
bayes_R2(m1)



# 5. Create plots ---------------------------------------------------------

# Use brms to compute dyadic social and partisanship effects
# conditioning on mean and reference values

me <- conditional_effects(m1,
                          effects = c("rel", "prt"))


# Create relationship plot

rel_plot <-
  me$rel %>% 
  mutate(
    effect1__ =
      effect1__ %>% 
      factor(
        levels =
          c(
            "Other",
            "Your spouse or partner",
            "Another relative",
            "A friend",
            "A neighbour or co-worker"
          ),
        labels = 
          c(
            "Other",
            "Your spouse\nor partner",
            "Another\nrelative",
            "A friend",
            "A neighbour\nor co-worker"
          )
      )
  ) %>% 
  ggplot(
    aes(
      x = estimate__,
      xmin = lower__,
      xmax = upper__,
      y = effect1__
    )
  ) +
  geom_errorbarh(height = 0,
                 size = .5) +
  geom_point(size = 1) +
  geom_point(size = .3,
             colour = "white") +
  geom_text(data = me$rel %>% 
              mutate(
                effect1__ =
                  effect1__ %>% 
                  factor(
                    levels =
                      c(
                        "Other",
                        "Your spouse or partner",
                        "Another relative",
                        "A friend",
                        "A neighbour or co-worker"
                      ),
                    labels = 
                      c(
                        "Other",
                        "Your spouse\nor partner",
                        "Another\nrelative",
                        "A friend",
                        "A neighbour\nor co-worker"
                      )
                  )
              ),
            aes(label = scales::percent(estimate__,
                                        accuracy = 0.1,
                                        suffix = ""),
                x = estimate__ + .002,
                y = effect1__),
            family = "Cabin",
            size = 2,
            vjust = 2,
            hjust = 0.7,
            inherit.aes = FALSE) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0, .35)) +
  labs(title = "Dyadic Relationship Type",
       y = "",
       x = "Prob. of Discussant Approval") +
  theme_bailey() +
  theme(legend.position = "bottom",
        legend.margin = margin(t = 0, b = 0, unit = "cm"),
        legend.title = element_blank(),
        legend.text = element_text(size = rel(1)),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(hjust= .5))
rel_plot


# Create partisanship plot

prt_plot <-
  me$prt %>% 
  mutate(
    effect1__ =
      effect1__ %>% 
      factor(
        levels =
          c(
            "Shared Non-partisan",
            "Shared Partisanship",
            "Opposing Partisanship",
            "Respondent Non-partisan",
            "Discussant Non-partisan"
          ),
        labels = 
          c(
            "Shared\nNon-partisan",
            "Shared\nPartisanship",
            "Opposing\nPartisanship",
            "Respondent\nNon-partisan",
            "Discussant\nNon-partisan"
          )
      )
  ) %>% 
  ggplot(
    aes(
      x = estimate__,
      xmin = lower__,
      xmax = upper__,
      y = effect1__
    )
  ) +
  geom_errorbarh(height = 0,
                 size = .5) +
  geom_point(size = 1) +
  geom_point(size = .3,
             colour = "white") +
  geom_text(data = me$prt %>% 
              mutate(
                effect1__ =
                  effect1__ %>% 
                  factor(
                    levels =
                      c(
                        "Shared Non-partisan",
                        "Shared Partisanship",
                        "Opposing Partisanship",
                        "Respondent Non-partisan",
                        "Discussant Non-partisan"
                      ),
                    labels = 
                      c(
                        "Shared\nNon-partisan",
                        "Shared\nPartisanship",
                        "Opposing\nPartisanship",
                        "Respondent\nNon-partisan",
                        "Discussant\nNon-partisan"
                      )
                  )
              ),
            aes(label = scales::percent(estimate__,
                                        accuracy = 0.1,
                                        suffix = ""),
                x = estimate__ + .002,
                y = effect1__),
            family = "Cabin",
            size = 2,
            vjust = 2,
            hjust = 0.7,
            inherit.aes = FALSE) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0, .35)) +
  labs(title = "Dyadic Party Identification",
       y = "",
       x = "Prob. of Discussant Approval") +
  theme_bailey() +
  theme(legend.position = "bottom",
        legend.margin = margin(t = 0, b = 0, unit = "cm"),
        legend.title = element_blank(),
        legend.text = element_text(size = rel(1)),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(hjust= .5))
prt_plot


# Group plots together

plot1 <- ggarrange(rel_plot, prt_plot, nrow = 1)


# Now we can render the plot as a jpeg file

jpeg(file = here("_output", "w12_approval.jpeg"),
     res = 300,
     width = 5.97,
     height = 2.2,
     units = "in")
plot1
dev.off()
